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\0 ' Abstract 



The orientational ordering transition is investigated in the quantum general- 
ization of the anisotropic-planar-rotor model in the low temperature regime. 



C. 

O , The phase diagram of the model is first analyzed within the mean-field ap- 

O 

proximation. This predicts at T = a phase transition from the ordered to 
the disordered state when the strength of quantum fluctuations, characterized 
by the rotational constant 0, exceeds a critical value . As a function of 
temperature, mean-field theory predicts a range of values of where the sys- 
tem develops long-range order upon cooling, but enters again into a disordered 
state at sufficiently low temperatures (reentrance) . The model is further stud- 
ied by means of path integral Monte Carlo simulations in combination with 
finite-size scaling techniques, concentrating on the region of parameter space 
where reentrance is predicted to occur. The phase diagram determined from 
the simulations does not seem to exhibit reentrant behavior; at intermediate 
temperatures a pronounced increase of short-range order is observed rather 
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than a genuine long-range order. 
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I. MOTIVATION 



Physisorbates are experimental realizations of quasi two-dimensional systems that dis- 
play an extremely rich phase behavior due to the competition between intermolecular and 
molecule-surface interactions, as documented, e.g., in Refs. [Q-Q]. Correspondingly there is a 
wealth of phase transitions between the various ordered phases as a function of temperature 
and coverage. Since many of these transitions occur at fairly low temperatures, quantum 
effects might play an important or even crucial role as recently demonstrated for the or- 
dering of hydrogen isotopes on graphite M . Molecular systems are particularly interesting as 
they possess orientational degrees of freedom which can order in addition to the positions || . 
In case of linear molecules, the anisotropic-planar-rotor (APR) model was devised to 
describe the herringbone (quadrupolar orientational two-sublattice) ordering transition ||, 
e.g., in commensurate N 2 monolayers on graphite. The classical two-dimensional APR model 
consists of planar rotators pinned with their center of rotation on a triangular lattice and 
interacting via nearest-neighbor quadrupolar interactions only; a three-dimensional version 
has also been proposed and investigated in various approximations ||. 

Over the years the APR Hamiltonian acquired the status of a statistical mechanical 
model on its own right, see Refs. flT0|,|5|j. Many of these activities arose because the order 



of the APR phase transition turned out to be extremely challenging to determine fi- 



nally favoring a first order phase transition that is "weak" and fluctuation driven [JL2 



sec 



the detailed review in Ref. f|. The plain APR model was generalized to include vacan- 
cies or impurities Hl3|-[T5[|, as well as multipole interactions other than those of quadrupolar 



symmetry |14,16fl. Simplified Hamiltonians were obtained after discretizing the continuous 



rotations to an anisotropic 6-state model |L7| or to more general discrete models |18fl . We 
introduced quantum generalizations of interacting two-dimensional quadrupolar lattice sys- 



tems |T9|j20|| , in particular of the APR model pT] , p2| . For a given lattice, the classical APR 
model has no free parameter because the quadrupolar coupling constant K only sets the 
energy scale of the model. Correspondingly, the phase diagram is one-dimensional and is 
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fully characterized by a single number, the transition temperature scaled by the constant K. 
In the quantum case, however, the quantum kinetic energy is determined by the mechanical 
moment of inertia I associated with the angular motion of the two-dimensional rotators. 
The resulting rotational constant G = h 2 /2I is an additional independent parameter that 
determines the strength and energy scale of the quantum fluctuations, and in the limit 
— > the quantum APR model reduces to the classical one. The resulting behavior of the 
quantum APR model is governed by the interplay between thermal and quantum fluctua- 
tions. This interesting feature is present also in other systems, such as, e.g., Ising models 
in transverse field |]23|j24l| , models for granular superconductors |25] and superconducting 



arrays p6 |, lattice 4 -theory p7j, and quantum four-state clock models | 28|| , see Ref. |22] for 



a short discussion of these related models. 



In our previous studies 1 2 1,22] we presented the first qualitative exploration of the 2D 



quantum APR Hamiltonian, using path integral Monte Carlo (PIMC) simulations [29 
adapted to rotational motion restricted to two dimensions []30| . It was demonstrated numer- 
ically [^2J that low-order approximation schemes such as quasiclassical Monte Carlo simula- 
tions using the quadratic Feynman-Hibbs effective potential, and the simple quasiharmonic 
approximation are useful only in the regime of small rotational constants G — > 0, whereas 
they fail completely in the large-0 range that is of interest here. The phase boundaries 
were estimated phenomenologically from the behavior of the order parameter, i.e., without 
applying any kind of finite-size scaling to the data obtained from 30 2 = 900 interacting 
rotators. Increasing the value of G, four distinct regimes were found based on the shape of 
the orientational order parameter. For small 0, the transition temperature and value of the 
ground state order parameter obtained from the classical model get only renormalized to 
smaller values. In the opposite limit of large B, the quantum fluctuations are that strong 
that they do not allow for any ordering even in the ground state. More interesting are the 
intermediate regimes. For increasing a crossover was found where there is residual ground 
state order, although significantly depressed from its classical value, but the order parameter 
inside the ordered phase first increases upon heating and goes through a maximum at some 



intermediate temperature before it decays in the disordered phase. For even larger 0, there 
seemed to be a region in parameter space with vanishing ground state order together with 
residual order at finite temperatures. The "tentative and qualitative phase diagram" in the 
Q-T plane thus seemed to exhibit a reentrance phenomenon in a range of where quan- 
tum and thermal fluctuations are competitive, see Fig. 9 in Ref. [2~2~|. A similar behavior 
was already discovered in mean- field studies of certain quantum Ising models PHI , quantum 
four-state clock models with quadrupolar interactions [28[] and interacting 3D quadrupolar 
rotators | 3lfl . The reentrance in the two-level systems was also found to be present when 
fluctuations were included to lowest order in form of a Kirkwood correction on top of the 



Bragg- Williams expressions |24j]. The reentrance in the three-dimensional system was as- 



cribed in Ref. |3T| to two distinct first order phase transitions, a standard order- disorder 
transition at higher temperatures that is also present in the corresponding classical model, 
and a transition of pure quantum nature at low temperatures. 

With the present paper, we concentrate on the region of the phase diagram of the 2D 
quantum APR model where reentrance might be expected to occur. We use first a mean- 
field approximation and then a PIMC simulation. The aim of the work is to study the 
orientational transition behavior in the reentrance regime. In particular, we are interested 
to find out whether there are actually two distinct phase transitions occurring as a function of 
temperature. To this end, we simulate the system on several length scales and use finite-size 



scaling, in particular Binder's cumulant |32],[33|, to analyze the order parameter distributions. 

The main body of the paper is organized as follows. The quantum APR model together 
with its order parameter is defined in Sect. ||. In Sect. [Ill A] , the quantum PIMC simulation 



technique is described, followed by the definition of the relevant observables in Sect. f \ 1 B 
and an outline of the finite-size scaling analysis method in Sect. |111 C. The mean-field 
approximation is worked out in Sect. [IV]. In Sect. [VJ, the results of the simulation are 
discussed and compared to those of the mean-field approximation and a possible scenario is 
suggested. In Sect. |VTJ, we summarize the results and draw some conclusions. 
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II. QUANTUM ANISOTROPIC-PLANAR-ROTOR MODEL 



Subject of this study are the properties and in particular the phase diagram of the 
quantum-mechanical version of the APR model The classical APR Hamiltonian was 

extensively studied over the years, both numerically and analytically, see Ref. || where 
the properties of the classical model including the order parameter and the herringbone 
orient at ionally ordered ground state are discussed. The N quadrupolar rotators are fixed 
with their center of mass on an ideal rigid triangular lattice {Rj} and only rotations in 
the two-dimensional surface plane are allowed. The interactions are truncated at the first- 
neighbor shell, and the rotators interact exclusively with their six nearest neighbors via 
the anisotropic part of a quadrupole-quadrupole potential of strength K in two dimensions. 
The quantum generalization results after supplementing this classical Hamiltonian with 
a non-commuting angular momentum part [Lj, (fi] = —ih5j t i which introduces quantum 
dispersion and thus qualitatively new effects due to additional fluctuations and tunneling. 
Correspondingly, the quantum APR Hamiltonian reads 

N r 2 N 

N Q2 N 

= - E ^-s + K E cos + 2 ^ - > ( 2 - 2 ) 

where the angle (pj of the j-th rotator pinned at site Rj is defined relative to one symmetry 
axes of the triangular lattice, and the six phases 0^, measure the angle between neighboring 
sites Rj and Rj on this lattice, i.e., 4>j,i {0, 7r/3, 27r/3, 7r, 47r/3, 57r/3}. We stress that 
for this system the triangular lattice structure is essential in order to produce a nontrivial 
ordered phase. On a simple square lattice the quadrupolar interaction would just favor 
unfrustrated perpendicular nearest-neighbor orientations of the characteristic T-shape. 

The moment of inertia I determines the rotational constant = h 2 /2I, which is the 
parameter that controls the strength of quantum effects. The other parameter of the model, 
which is the quadrupolar coupling constant K, can be conveniently taken as the energy 
and temperature scale. We can thus reduce all quantities related to energies by K, and 
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define, e.g., the dimensionless temperature T* = ksT/K, energy E* = E/K, and rotational 
constant 0* = Q/K. 

The long range order parameter sensitive to herringbone ordering of the rotational axis 
has three (a = 1, 2, 3) independent components 

1 N 

$ ° = J7 E sin ( 2 <^ -2r/ a )exp[iQ Q R i ] , (2.3) 

3=1 

where 

Qi = tt(0,2/V3), vi = 0, 

Q 2 = 7r(-l,-l/V3), m = 2tt/3, 

Q 3 = tt(1,-1/V3), r/ 3 = 4vr/3. (2.4) 



III. NUMERICAL METHODS 



A. Path Integral Monte Carlo Simulation Technique 



We study the properties of the Hamiltonian ( |2.2j ) by means of PIMC simulations, for 
general reviews on the concept of PIMC we refer to Ref. P5| . We just stress here that for 
rotational motion in two dimensions new features connected to the restricted integration 
space show up in the formalism [33]. This latter aspect and our implementation of these 
specialties in a PIMC scheme is discussed in detail in Ref. [jn|, and here we present only 
the essential features. 

In general the partition function of a Hamiltonian of type (|2.1|) is given by 



Z = Tr exp \-(3H] 



( IP 
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where V{<Pj\ iff ) denotes the APR pair potential of (|2.2|) evaluated separately for the con- 



3 

figuration at each imaginary time slice s = 1, • • • , P. Each quantum mechanical rotational 
degree of freedom is represented in this path integral representation by P classical rotators 
which form closed loops and interact via harmonic type interactions; for the related ring- 



polymer picture see Ref. [35,29]. The parameter P is the Trotter number, the configuration 
{(pf \ (ff, ■ • - , <Pj } is a realization of a Trotter path, and the path integral results from 
the proper integration and summation over all possible paths. However, contrary to path 
integrals for translational degrees of freedom these loops need not to be closed using periodic 
boundary conditions, but only modulo 2ir; note that the classical angles tpf of (|3.2|) are 
not confined to [0, 2ir) but are allowed on the whole interval [—00,00]. The resulting mis- 
match rij is called the "winding number "of the j th path |34| and the formulation (|3.2|) is the 
"winding number representation" of the partition function. Only the Boltzmann- weighted 
summation over all possible winding numbers in addition to the integration over all paths 
having a certain winding number yields the correct quantum partition function in the Trot- 



ter limit P — > 00; see Ref. JT9| for a full discussion of that issue. Thus we have to include in 
the algorithm in addition to local and global moves of the angular degrees of freedom {<£>j- } 
also attempts to change the winding numbers {rij} of the individual rotators. Our algorithm 



was tested p3| , |T9"l| against exact (single-particle) results and a close agreement was observed. 
In a previous study |22| it turned out that the use of a finite Trotter number of P = 500 at a 
temperature of T* = 0.03 and a progressive linear decrease of P with increasing temperature 
was sufficient to be in the large P limit as required in ( |3.2|) . 

In the present paper we study in particular finite-size effects on the quantum APR 
transition. The required different linear dimensions L of the system were chosen to range 
from N = L 2 = 144 to 900, where additional data from Ref. J22| have been included in the 
analysis. We equilibrated the systems starting from the ideally ordered herringbone ground 
state and made where possible use of previous runs either at lower temperatures or smaller 
rotational constant. Our statistics is based on the order of 10 5 PIMC steps for each data 



S 



point, a typical CPU time per data point was of the order of 50 hours on a CRAY-YMP. 



B. Observables and Estimators 

The quantities we determined from the simulations are kinetic and potential energy, 
order parameter, quantum librational amplitudes, as well as both fourth order p. 11]) and 
second order ( |3.12| ) cumulants. 

The energies were obtained from the primitive estimator which proved to be sufficient 
in the present case |B|. The kinetic energy E kin and the potential energy E pot per particle 



£m» = ^ - (± £ £ 4^f - ¥>r " + ^„ P f) , (3.3) 



are given by 

Pk B T _ 

^ P ot = (^E^ E n^ s) )) , (3-4) 

s — 1 <j,i> 

and the total energy £tot is given by the sum of -Ekin and -Epot- The estimator for the order 
parameter component <3> Q is given by the expression 

$ - = W5 E E sin(2^ s) - 2 Va ) expfiQ^R,-] , (3.5) 

which trivially follows from ( |2.3| ). In the following we will use the total long range order 
parameter defined as the length []T§| , [I2" | 



* = (3-6) 



a=l 



of the three component vector order parameter (3T5). 



A quantity which measures the quantum delocalization of the rotational degrees of free- 



dom can be defined by the expression [JL9 



Xiv j=l r s=l r k=l 1 

This average spread is by its definition zero for a classical system and is thus a measure of 
the pure quantum contribution to the librations or rotations of the individual rotators. 
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C. Finite-Size Scaling 



A fingerprint of a continuous phase transition is the divergence of the correlation length 
£ at the critical temperature T c , with £ ~ \t\~ u for t := 1 — T/T c , \t\ <C 1. Under these 
conditions on the reduced distance from the phase transition t > the order parameter $ 
in the infinitely large system depends on t as $ ~ {t^, where v and /? denote here the usual 
critical exponents. In finite systems of linear dimension L the order parameter is given by 
the expression |B3| 



(<% = , (3.8) 

where $ is the scaling function that is associated to $. The order parameter distribution 
function Pl(&) in the finite system is thus a function of the scaling variables LP/ V Q and 

LIZ, 

P L (<$>) = L / v P{LP/ v $,L/£) , (3.9) 

and the fe-th moment of the order parameter distribution function is given by 

= lfl v J d$ $ fe P(L^$, L/i) = L- ph l v $ k (L/g) , (3.10) 

where the symbols ( ) L denote canonical averages with systems of linear dimension L. 

A very useful quantity for the determination of a critical point which is directly based 
on order parameter moments is the fourth order cumulant [|3^,[3t| or the second order 
cumulant |Q defined as 

U?> = 1- ^\ , (3.11) 

U® = l-^\, (3.12) 

where one can see that the explicit dependence on system size drops out if the k-th moments 
~ L~P k ' v are re-expressed with the aid of (|3.10|) in terms of their scaling functions. We only 
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compile here the main features and usage of these quantities, and refer for further details to 
the literature p2| , |33 |. 



In the case of second order transitions, the cumulants adopt a non-trivial universal value 
U* at the critical point, irrespective of system sizes {L} in the scaling limit. Thus, plotting 
Uj^\t) or Ul\t) for different linear dimensions as a function of temperature yields an 
intersection point U(T C ) =: U* which gives an accurate estimate of the critical temperature 
in the infinite system for a temperature driven second order transition. Below and above 
the transition, the cumulants flow to trivial limiting values depending on the definition 
of the order parameter; the larger the system, the faster the convergence. Instead of the 
temperature, other parameterizations of approaching a critical point can also be chosen. 
More recently, it was worked out that the concept of order parameter cumulant crossings 



is also useful to analyze first order phase transitions ||37|| . In this case, one can observe an 
effective crossing point at a non-universal value U* at the phase transition. The approach 
of both the transition point and the value of U at the transition to the infinite system 
limit is quite fast: the correction depends roughly speaking on the inverse volume of the 
system J57[] . Thus, for practical numerical purposes the order parameter cumulant can be 
taken as acquiring an intersection point at a first order transition in a way similar to that 
occurring at a second order transition. 

IV. MEAN-FIELD THEORY 



In this section we start from the Hamiltonian ( |2.2|) and determine the phase diagram for 
our model using a mean-field approximation. This consists in considering a single quantum 
rotator in the mean field of its six nearest neighbors and finding a self-consistent condition 
for the order parameter. Solving the latter condition, we determine the phase boundary and 
also the order of the transition. Our mean-field approximation is similar in spirit to that 



used in Ref. 31] for the case of 3D rotators. 



We shall assume that the order parameter component $1 becomes non-zero in the ordered 
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phase, $1 = ±(sin2</?) 7^ (we choose positive sign for the single rotator, note that the sign 
alternates passing from one row of the herringbone structure to another, see Eq. fl2.3| )), while 
$2 = $3 = 0. Now it is useful for a moment to come back to normal (non-reduced) units for 
G, K and T. In order to proceed, we write each interaction term K cos(2<£> + 2cp' — 40^) as 
a product of trigonometric functions depending separately on the angular variables p of the 
single rotator and p> of its nearest neighbors. Averaging over the variables p> we find the 
following contributions: for the two terms corresponding to <f>ij = we get — 2K<&\ sin 2p>, 
for the two terms corresponding to <fiij = n/3 we get — KQ\ sin2y? + \^3K&i cos2y) and for 
the two terms corresponding to <fiij = 27r/3 we get — K$>i sin 

2y> - v^-RT^i cos 

2<p. Summing 

these contributions from all six nearest neighbors, we find the total mean-field potential 
acting on the single rotator, which reads 

H%* = -4K* 1 sin 2<p . (4.1) 

Adding the kinetic energy, we get the corresponding one-particle Schrodinger equation for 
the on-site problem 

H MF ^ = ( -0 ^— - 4^$! sin2y? ) ^ = . (4.2) 
\ dtp 2 J 

We now introduce the quantity q = —2K<&i/Q and perform a trivial shift of the angular 
variable <p by defining 9 = p — ir/4. In terms of these new variables, the eigenvalue problem 
04.2| ) can be written as 

^ + (|-2gcos2^)vl/ = 0, (4.3) 

which is the well-known Mathieu's equation ||38|| . Its eigenvalues can be labeled by an non- 
negative integer number m, and for each m^O there are two eigenvalues, a m (q) associated 
with an even periodic solution, and b m (q) associated with an odd periodic solution. For 
m = there is just one eigenvalue ao(q) associated with an even periodic solution. 

In order to find the self-consistent condition, we have to determine the order parameter 

= (4 ,) 
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as a function of the parameter q. In the last equation we have used the free energy per 
site defined as f = —{1/(3) lnX^e - ^. Solving the equation ( |4.4| ) simultaneously with the 
condition q = —2K$i/Q, we find the equilibrium value of the order parameter for given 
values of inverse temperature /3 and model parameters K, 0. 

Because our main interest is to determine the phase boundary and the order of the 
transition, we do not have to find the complete expression for <&i(q). For a continuous phase 
transition, the phase boundary is a curve in parameter space on which a nontrivial and 
infinitesimally small solution $1 ^ appears, apart from the trivial one $i = 0, which is 
always present and corresponds to the disordered phase. In order to determine the phase 
boundary and check the order of the transition, it is enough to find an expansion of $i in 
powers of q up to the third order, which can be obtained from the expansion of fo in powers 
of q up to the fourth order. The expansions of the eigenvalues of the Mathieu equation can 



be found in [35] and we quote here just first few of them, which have the following form 



q 2 7o 4 

a »<'> = -¥ + lf8-- < 4 ' 5 > 

°i(-*) = M*) = l-?- ¥ + 64 - 1536-- • 

For practical calculations in the low temperature region, we can truncate the infinite sum 
in /o and include only a finite number of terms. Increasing the number of terms would yield 
a progressively better approximation in the high temperature limit, but this is actually not 
necessary, since in this limit the quantum effects become negligible and the model behaves 
classically. The classical model has already been investigated within various mean-field 



approximations in Ref. [p9|| . We have taken the eigenvalues corresponding torn < 6. The 
first neglected eigenvalue then equals 7 2 = 49 for q = and therefore our approximation 
should yield reliable results up to temperatures which are at least an order of magnitude 
lower than the first neglected term, therefore for temperatures T/Q < 5. Substituting 
the energy eigenvalues expansions into ( |4.4| ) and calculating the free energy fo and order 
parameter <3>i (the actual symbolic calculation has been performed by Mathematica) , we find 
that the expansion of the latter has the form 
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$1 = Xi<? + X3<? 3 + ••• , (4.6) 

where the coefficients xi an d Xs are somewhat complicated expressions. Combining the 
expression (|4.6[ ) with the equation g = — 2A"$i/0, we find that the condition for the phase 
boundary is given by K c = — 0/(2%i). We quote here the final result for K c , which reads 

Kc = 420 *< Z * (47) 

420 + 2105! + l) - 280fi 2 - 105^3 - 565 4 - 35B 5 - 245 6 ' 

where B m = exp[-m 2 c /T c ] and Z 6 = Em=-e B m- 

We yet have to make sure that the transition is really continuous. The necessary and 
sufficient condition for this is that the coefficient X3 is negative on the phase boundary ( [4.71 ). 
We have found this to be always satisfied. Therefore the phase boundary Q4.7| ) corresponds 
to a continuous phase transition. In the high-temperature, classical limit this agrees with 
the finding in Ref. f39 |. In order to map the phase diagram in the 0*-T* plane, we have to 



set K c = 1 in (|4.7| ) and solve for T* as a function of 0*. This can be done numerically and 
the resulting phase diagram is shown on Fig. [|. 

V. RESULTS AND DISCUSSION 

To start the discussion, we note that the present study of the quantum APR model was 
partly motivated by the strong changes in shape of the orientational order parameter <3> as 
a function of temperature as the rotational constant was increasing from its classical value 



0* = 0, see Fig. 3 in Ref. p2| . For small enough 0* it was found that the order parameter 
decays monotonously with increasing temperature similarly to the classical case. This is 
qualitatively different for larger 0*, where $(T*) becomes a non-monotonous function of 
temperature. In the present paper, we study this regime in much more detail and vary in 
addition to 0* also the linear dimension L of the lattice. 

In Fig. [I], the order parameter vs. temperature curves are shown for different system 
sizes and for two representative 0* choices, 0* = 0.6109 in (a) and 0.6982 in (b). As can 
be seen from the data, in case (a) the order parameter at low temperatures has a non-zero 
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value and on heating first increases up to a maximum and then decreases, which means 
that the system develops most order at intermediate temperatures. The results obtained for 
different system sizes indicate that in the high temperature region the order parameter scales 
to zero with increasing system size, whereas the maximum value is system size independent. 
Therefore there is a long-range ordered phase at low temperatures present in addition to 
the disordered high temperature phase. This orientational ordering that takes place with 
decreasing temperature results in the limit 0* — > in the classical APR transition and 
the well-defined herringbone ordered ground state. This picture changes, however, as the 
quantum fluctuations get more pronounced. When the value of Q* is increased further, 
the maximum value of the long-range order parameter drops considerably with increasing 
system size, as can be seen for G* = 0.6982 from Fig. |] (b). The order parameter at our 
lowest temperatures is now strongly depressed and also decreases further with increasing 
linear dimension of the system. The finite-size effects in this intermediate Q* range are thus 
pronounced not only at the high-temperature wing of the order parameter vs. temperature 
curve, but at all temperatures and in particular also close to the maximum of the order 
parameter. This is already a hint that the residual order in this range of G*-values might 
actually be a finite-size artifact. From the previous study |2!| we know that a further 



increase in G* would suppress the ordering even more and finally lead to a situation with 
vanishing order parameter at any temperature. All this motivated us to analyze the behavior 
of the system further in order to see whether the intermediate-temperature maximum and 
low-temperature reduction of ordering, suggestive of a reentrant orientational melting, are 
really associated with the latter phenomenon. In order to answer such question one has 
to employ finite-size scaling techniques which would allow to find out whether a particular 
phase transition is present, and if so, to locate it. 

Before we start with the numerical finite-size scaling analysis of the PIMC data, let us 
pause for a moment and discuss the mean-field predictions from the previous section. The 
resulting phase diagram in the G*-T* plane is shown in Fig. |2| with the solid line. We note 
that in the classical limiting case G* — > we find T* = 2 which agrees with earlier mean- 
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field calculations |39| . Any crossing of the phase boundary in the mean-field phase diagram 



Fig. ^] corresponds to a second order phase transition; note that it is believed that for 

O* = 0, the classical APR model undergoes a fluctuation driven weak first order transition 
at T* « 0.76 for L — > oo. We see that at zero temperature there is a quantum phase tran- 
sition at the value of Q* MF = 1. The most interesting feature of the phase diagram is that 
there is a region of rotational constants ranging from 1 to roughly 1.25 for which the system 
is ordered at intermediate temperatures but disordered in the ground state (reentrance) . 
The intuitive explanation of this phenomenon is the following. At low temperatures, the 
individual rotors are mostly in their totally rotationally symmetric ground state, which does 
not possess quadrupolar moment and therefore cannot induce ordering via the quadrupolar 
term. At intermediate temperatures, the excited states with non-zero quadrupolar moment 
become populated and induce ordering which persists to larger values of the rotational con- 
stant. According to the mean-field theory, reentrance takes place for a rather broad range 
of rotational constants 1 < 0* < 1.25. This corresponds to a decrease of the critical ro- 
tational constant by roughly 20% from its maximum value of about 1.25 to the value of 1 
at the ground state transition, and represents a well pronounced feature. Concerning the 
validity of this mean-field result, it is known that the mean-field approximation sometimes 



tends to overemphasize or even create reentrant behavior, as pointed out in Ref. [p4j| . On 
the other hand, reentrant behavior has been experimentally observed in solid HD under 
compression f40|| . This three-dimensional system, however, although consisting of diatomic 
molecules interacting via approximate quadrupolar interactions, differs fundamentally from 
our model in spatial dimensionality and structure of the lattice as well as in the dimension- 
ality of the order parameter. In order to settle the question of reentrance in the 2D quantum 
APR model, we present now the results of numerical simulations and analysis techniques. 

We start with the results for potential, kinetic and total energies of our model. These 
are presented in Figs. [| and f| for two representative cases, 9* = 0.6109 and 0* = 0.6982, 
respectively. In both cases the potential energy first decreases with temperature until maxi- 
mum order is achieved, while the kinetic energy increases strongly in this region. For larger 
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temperatures both E pot and E^ iD increase with temperature. Thus, when the temperature 
is increased from zero the rotators occupy higher rotational states which allow for more 
pronounced orientational ordering compared to that in the ground state. In these ordered 
states the attractive quadrupolar interaction is larger resulting in a lower potential energy, 
and the kinetic energy increases strongly with temperature due to enhanced localization of 
the rotators along some direction. This scenario continues with temperature until maximum 
order is achieved. Further increase of temperature results in an increase of E pot due to ther- 
mal disorder making the quadrupolar interaction less important compared to the thermal 
energy. This can finally lead to a phase transition from the low temperature orient at ionally 
ordered phase to a high temperature disordered phase, provided the rotational constant 
does not exceed a certain value. At very low temperatures the slope of the total energy as 
a function of temperature, i.e., the specific heat, approaches zero, as expected for quantum 
systems. We did not find a strong size dependency of E tot (T*) and thus the specific heat 
behavior also does not seem to depend on the system size. 

Similar to the total energy, the average spread R^,, which is a measure of the quantum 
mechanical derealization of the rotators, is not dependent on system size as shown in Fig. 
As can be seen from its definition ( |3.7| ) this quantity is a single particle property that is 
by construction not particularly sensitive to collective effects. The spread approaches its 
classical limit, i.e., zero, for large temperatures, the approach being slower for larger 0*. 
In the limit of low temperatures, it reaches a ground state value that is larger than 90° 
for the rotational constant G* = 0.6982, which means that an individual rotator is no more 
confined to perform librational motion around a preferred orientation but is instead strongly 
delocalized. On the other hand we know already from Fig. [I] (b) that the order parameter is 
vanishingly small at low temperatures and decreases even more with increasing system size 
for this value of 0*. The behavior of the average spread is thus a clear demonstration that it 
is the quantum tunneling which induces the disordering of the ground state for sufficiently 
large rotational constants. 

In order to finally address the question whether our system has a reentrant phase tran- 
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sition as predicted by the mean-field study we analyzed the low temperature region by the 
cumulant intersection finite-size scaling method described in Sect. |111 C| , see Figs. |6], [7], and 



|8]. For our smallest rotational constant 0* = 0.6109 we clearly find an intersection point for 
both U ( l ] and uf ] at about T* « 0.303, see Figs. |6](a) and 0(a). For the larger rotational 
constant G* = 0.6364 (b), the intersections of the fourth and second order cumulants occur 
again both at the same value within the error bars. These crossings arise because the cumu- 
lants for larger systems approach their limiting values faster than those for smaller systems. 
For Q* = 0.6666 (c) the cumulants on the different length scales at low temperatures have 
values which cannot be distinguished within the error bars, see Figs. ||(c) and 0(c). The 
large value of U* ~ 0.65 that was found for the classical APR model |12| causes here the 



problem that within our numerical accuracy we cannot identify an intersection point but 
rather obtain a whole temperature region that is characterized by pronounced fluctuations. 
For larger G* constants the behavior of the two cumulants can again be distinguished fairly 
well. However, contrary to what we found for G* < 0.6364, the larger system has now the 
smaller cumulant throughout the whole temperature range, see Figs. ^|(d) and 0(d). This 
signals the presence of the orient at ionally disordered phase for G* = 0.6982 that extends 
from the ground state up to high temperatures. From this behavior of the cumulants we 
conclude that there is no evidence for the reentrance transition but rather a pronounced 
increase of short-range order at intermediate temperatures in the neighborhood of the phase 
boundary line. The latter has at low temperatures a roughly vertical (i.e. Q*-independent) 
slope, rather than the characteristic shape suggested by the mean-field result in Fig. |2|. 

In order to reinforce such conclusion for this part of the phase boundary we also studied 
the cumulants as functions of the rotational constant at constant temperature in the range 
T* < 0.2424. The fourth order cumulant plots of Fig. |8| show that the intersection occurs 
in this low temperature range at a value of about Q* ~ 0.667 ± 0.015, which stays constant 

(2) 

within the numerical noise; the second order cumulants U L show the same behavior. Note 
that a systematic increase followed by a decrease of this quantity is expected as a function 
of temperature if reentrance does occur. We can furthermore infer that the nontrivial value 



U* of the cumulants at crossing is also in this low temperature regime with pronounced 
tunneling and quantum effects within the numerical accuracy identical to the value U* ~ 0.65 
found for the classical APR phase transition ||12|| . We are thus lead to conclude that the 
APR transition temperature decreases slowly from its classical limit value of T* « 0.76 at 
0* = down to about T* « 0.24 at O* ~ 0.67, where it suddenly drops dramatically. 
This numerically obtained non-reentrant phase diagram is included in Fig. |^. The PIMC 
simulation results clearly exclude the existence of a strongly pronounced reentrant feature 
in the phase diagram. Of course, we cannot exclude the possible existence of a narrow 
reentrance region falling within the error bars of the present data. However, in any case 
these error bars are much smaller than the 20 % decrease of the critical rotational constant 
predicted by the mean-field theory. The latter approximation is successful in predicting the 
phenomenon of enhanced ordering at intermediate temperatures, but the deficiency is that 
it exaggerates the range of order and incorrectly predicts it to become long-ranged. As can 
be seen from Fig. |2|, the mean-field theory apparently treats the quantum fluctuations (the 
limit T* — > 0) better than the thermal fluctuations (the limit 0* — > 0), since the transition 
point of the quantum induced transition is overestimated only by a factor of about 1.5, 
whereas the purely classical transition is off by a factor of more than 2. 



VI. SUMMARY AND CONCLUSIONS 

We have studied the quantum generalization of the anisotropic-planar-rotor model con- 
sisting of point quadrupoles that are pinned with their center of rotation on a triangular 
lattice. This two-dimensional statistical mechanical model system exhibits interesting ori- 
entational ordering effects, as a function of temperature and the rotational constant which 
controls the strength of the quantum fluctuations. A mean-field study predicts for this quan- 
tum system a reentrant orientational order- disorder transition with a pronounced reentrance 
region. Thus, in an intermediate range of rotational constant values at low temperatures 
the system orders with increasing temperature, and then disorders again at higher tempera- 



19 



tures via a phase transition, which corresponds to the well-known herringbone orientational 
transition in the classical limit. 

We studied the quantum APR model also numerically using path integral Monte Carlo 
simulations in combination with a finite-size scaling analysis in order to explore the region of 
the phase diagram where reentrance is expected. It turned out that no reentrant transition 
is present, at least within the usual numerical limitations of such simulations, which are 
most importantly set by the statistical errors and limited system sizes concerning both the 
N and P dimensionalities. In order to determine whether there still is a tiny reentrance 
region hidden within our numerical accuracy, at least an order of magnitude larger amount 
of computer time would be required, which at present seems to be prohibitively expensive. 
For the same reason we did not attempt to address the issue of the order of the quantum 
induced APR phase transition in the ground state. We could, however, infer that the non- 
trivial cumulant values at the APR phase transition at low temperatures and large rotational 
constants are within the estimated uncertainty identical to the value obtained in the classical 
limit. 

There are several possibilities to extend the present study. Concerning the problem of 
reentrance in systems consisting of quantum rotators, it might be interesting to study the 
effect of higher multipole interactions. In addition, it would be desirable to go beyond the 
mean-field approximation by means of analytical techniques, such as, e.g., renormalization 
techniques capable of including accurately the quantum fluctuations. Finally, it is certainly 
desirable to develop even more efficient quantum simulation and analysis techniques. All this 
demonstrates that investigations of phase transitions in quantum systems with continuous 
degrees of freedom are even today still a challenge. 
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FIGURES 

FIG. 1. Orientational order parameter $ as a function of temperature. Symbols indicate 
different system sizes, statistical error bars are shown, lines are for visual help only. The values of 
the rotational constants are: (a) 0* = 0.6109, (b) 6* = 0.6982. 

FIG. 2. Phase diagram of the quantum APR model in the 0*-T* plane. The solid curve 
shows the line of continuous phase transitions from an ordered phase at low temperatures and 
small rotational constants to a disordered phase according to the mean-field approximation. The 
symbols show the transitions found by the finite-size scaling analysis of the path integral Monte 
Carlo data. The dashed line connecting these data is for visual help only. 

FIG. 3. Energies as a function of temperature for the rotational constant 0* = 0.6109. (a) 
potential energy E pot , (b) kinetic energy E^ in , (c) total energy E tot . Symbols indicate different 
system sizes, lines are for visual help only. 

FIG. 4. Energies as a function of temperature for the rotational constant 0* = 0.6982. (a) 
potential energy E po t, (b) kinetic energy E^ n , (c) total energy .Etot- Symbols indicate different 
system sizes, lines are for visual help only. 

FIG. 5. Average spread R v in degrees as a function of temperature. The values of the 
rotational constants are: (a) 0* = 0.6109, (b) 0* = 0.6982. Symbols indicate different system 
sizes, lines are for visual help only. 

FIG. 6. Fourth order cumulant as a function of temperature. The values of the rotational 
constants are: (a) 0* = 0.6109, (b) 0* = 0.6364, (c) 0* = 0.6667, (d) 0* = 0.6982. Symbols 
indicate different system sizes, lines are for visual help only. 

(2) 

FIG. 7. Second order cumulant as a function of temperature. The values of the rotational 
constants are: (a) 0* = 0.6109, (b) 0* = 0.6364, (c) 0* = 0.6667, (d) 0* = 0.6982. Symbols 
indicate different system sizes, lines are for visual help only. 
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FIG. 8. Fourth order cumulant as a function of rotational constant. The values of the 
temperature are: (a) T* = 0.1515. (b) T* = 0.1818, (c) T* = 0.2121, (d) T* = 0.2424, Symbols 
indicate different system sizes, lines are for visual help only. 
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